/*To calculate first-order gradient in adjoint method
*
	Chunli Zhang, clzhang@mail.ustc.edu.cn, zwlab
    Last revised on July 27th, 2018*/
#include<stdio.h>
#include<math.h>
/* 
*Input:
  */
//2D case
int first_order_gradiant_ela(double *vx_inc, double *vz_inc, 
				double *tauxx_inc, double *tauzz_inc, double *tauxz_inc, 
				double *vx_adj, double *vz_adj, 
				double *tauxx_adj, double *tauzz_adj, double *tauxz_adj,
				double *vp, double *vs, double *rho, 
				double *grad_vp, double *grad_vs, double *grad_rho,
				int Xn, int Zn, int Tn)
{
	double *vx_corr, *vz_corr, *tauxx_corr, *tauzz_corr, *tauxz_corr;
	int num_grd= Xn*Zn;
	vx_corr=(double*)malloc(num_grd*sizeof(double));
	vz_corr=(double*)malloc(num_grd*sizeof(double));
	tauxx_corr=(double*)malloc(num_grd*sizeof(double));
	tauzz_corr=(double*)malloc(num_grd*sizeof(double));
	tauxz_corr=(double*)malloc(num_grd*sizeof(double));

	/*check the input parameters*/
	/*check the file*/
	/*calculate the correlations*/
	ela_correlation(double *vx_inc, double *vz_inc, double *tauxx_inc, double *tauzz_inc, double *tauxz_inc, 
				double *vx_adj, double *vz_adj, double *tauxx_adj, double *tauzz_adj, double *tauxz_adj, 
				double *vx_corr, double *vz_corr, double *tauxx_corr, double* tauzz_corr, double* tauxz_corr,
				int Xn, int Zn, int Tn)
	/*calculate gradient*/
	ela_gradient(double*vs,double*vp,double*rho,double*vx_corr,double*vz_corr,
				double*tauxx_corr,double*tauzz_corr,double*tauxz_corr,
				double*grad_vp,double*grad_vs, double*grad_rho, 
				int num_grd, double dt) 
	return 0;	    	
}
int first_order_gradiant_acu(double *vx_inc, double *vz_inc, 
				double *pressure_inc, 
				double *vx_adj, double *vz_adj, 
				double *pressure_adj,
				double *vp, double *rho, 
				double *grad_vp, double *grad_rho,
				int Xn, int Zn, int Tn)
{
	double *vx_corr, *vz_corr, *tauxx_corr;
	int num_grd= Xn*Zn;
	vx_corr=(double*)malloc(num_grd*sizeof(double));
	vz_corr=(double*)malloc(num_grd*sizeof(double));
	pressure_corr=(double*)malloc(num_grd*sizeof(double));
	/*check the input parameters*/
	/*check the file*/
	/*calculate the correlations*/
	ac_correlation(double *vx_inc, double *vz_inc, double *pressure_inc, 
				double *vx_adj, double *vz_adj, double *pressure_adj,  
				double *vx_corr, double *vz_corr, double *pressure_corr, 
				int Xn, int Zn, int Tn)
	/*calculate gradient*/
	ela_gradient(double *vp,double *rho,
				double *pressure_corr,double *vx_corr, double *vz_corr,
				double *grad_vp,  double *grad_rho, 
				int num_grd, double dt) 
	return 0;	    	
}

/*Function: correlation
 * input: 
 */
void ela_correlation(double *vx_inc, double *vz_inc, double *tauxx_inc, double *tauzz_inc, double *tauxz_inc, 
				double *vx_adj, double *vz_adj, double *tauxx_adj, double *tauzz_adj, double *tauxz_adj, 
				double *vx_corr, double *vz_corr, double *tauxx_corr, double* tauzz_corr, double* tauxz_corr,
				int Xn, int Zn, int Tn)
{
	int i,j;
	int num_grid;
	num_grd=Xn*Zn;
	corr(vx_inc,vx_adj,num_grd,Tn,vx_corr);
	corr(vz_inc,vz_adj,num_grd,Tn,vz_corr);
	corr2(tauxx_inc,tauxx_adj,tauzz_inc,tauzz_adj,num_grd,Tn,tauxx_corr,tauzz_corr);
	corr(tauxz_inc,tauxz_adj,num_grd,Tn,tauxz_corr);
	return;
}
void ac_correlation(double *vx_inc, double *vz_inc, double *pressure_inc, 
				double *vx_adj, double *vz_adj, double *pressure_adj, 
				double *vx_corr, double *vz_corr, double *pressure_corr, int Xn, int Zn, int Tn)
{
	int i,j;
	int num_grd;
	num_grd=Xn*Zn;
	corr(vx_inc,vx_adj,num_grd,Tn,vx_corr);
	corr(vz_inc,vz_adj,num_grd,Tn,vz_corr);
	corr(pressure_inc,pressure_adj,num_grd,Tn,pressure_corr);
	return;
}
/*Function: corr
 * cross correlation with zero lag */
void corr(double*array1, double*array2, int num_grd, int timestep, double*xcorr)
{
	int i,j;
	for(j=0;j<num_grd;j++)
	{
		xcorr[j]=0.0;
		for(i=0;i<timestep;i++)
		{
			xcorr[j]=xcorr[j]+array1[i*num_grd+j]*array2[i*num_grd+j];
		}
	}
	
}
void corr2(double*array1, double*array2, double*array3, double*array4, int num_grd, int timestep, double*xcorr1, double*xcorr2)
{
	int i,j;
	for(j=0;j<num_grd;j++)
	{
		xcorr1[j]=0.0;
		xcorr2[j]=0.0;
		for(i=0;i<timestep;i++)
		{
			xcorr1[j]=xcorr[j]+(array1[i*num_grd+j]+array3[i*num_grd+j])*(array2[i*num_grd+j]+array4[i*num_grd+j]);
			xcorr2[j]=xcorr[j]+(array1[i*num_grd+j]-array3[i*num_grd+j])*(array2[i*num_grd+j]-array4[i*num_grd+j]);
		}
	}
	
}
/*Function: mem_alloc:Memory location (Temporary)
 * input:
 * */
void ela_gradient(double*vs,double*vp,double*rho,double*vx_corr,double*vz_corr,
				double*tauxx_corr,double*tauzz_corr,double*tauxz_corr,
				double*grad_vp,double*grad_vs, double*grad_rho, 
				int num_grd, double dt) 
{
//
	int i,j;
	double *lammu2, *mu2, *grad_mu, *grad_lam, *grad_rhos;
	//*allocate local parameters
	lammu2=(double*)malloc(num_grd*sizeof(double));
	mu2=(double*)malloc(num_grd*sizeof(double));
	grad_mu=(double*)malloc(num_grd*sizeof(double));
	grad_lam=(double*)malloc(num_grd*sizeof(double));
	grad_rhos=(double*)malloc(num_grd*sizeof(double));
	
	for(j=0;j<num_grd;j++)
	{
		lammu2[j]=pow((lam[j]+mu[j]),2);
		mu2[j]=pow(mu[j],2);
		grad_lam[j]=-dt*(0.25*tauxx_corr[j]/lammu2[j]);
	}
	for(j=0;j<num_grd;j++)
	{
    	grad_mu[j]=-dt*(0.25*tauxx_corr[j]/lammu2[j]+0.25*
			tauzz_corr[j]/mu2[j]+tauxz_corr[j]/mu2[j]);
		grad_rhos[j]=-dt*(vx_corr[j]+vy_corr[j]);
		grad_vp[j]=2.0*vp[j]*rho[j]*grad_lam[j];
		grad_vs[j]=-4.0*vp[j]*rho[j]*grad_lam[j]+2.0*rho[j]*vs[j]*grad_mu[j];
		grad_rho[j]=(pow(vp[j],2)-2.0*pow(vs[j],2))*grad_lam[j]+pow(vs[j],2)*grad_mu[j] +grad_rhos[j];
	}
}
void acu_gradient(double*vp, double*rho, 
				double*pressure_corr, double*vx_corr, double*vz_corr,
				double*grad_vp, double*grad_rho, int num_grd, double dt)
{
	int i,j;
	double *lam2, *grad_rhos, *grad_lam;
	lam2=(double*)malloc(num_grd*sizeof(double));
	grad_rhos=(double*)malloc(num_grd*sizeof(double));
	for(j=0;j<num_grd;j++)
	{
		lam2[j]=pow(lam[j],2);
		grad_lam[j]=-dt*(0.25*pressure_corr[j]/lam2[j]);
		grad_rhos[j]=-dt*(vx_corr[j]+vz_corr[j]);
	}
	for(j=0;j<num_grd;j++)
	{
		grad_vp[j]=2.0*vp[j]*rho[j]*grad_lam[j];
		grad_rho[j]=(vp**2)*grad_lam[j] + grad_rhos[j];
	}
}
/*Function: mem_alloc:Memory location (Temporary)
 * input:
 * */
